Twittter Github Email

Course materials for 2020-11-2 AFEC at XTBG.

1 Prerequisites

install.packages("spdep")
library(tidyverse)
library(spdep)

2 Spatial autocorrelaiton

  • Neighbors affect the observation.

  • Let’s condier maps with some environmental variables (e.g., soil N, soil moisture…).

  • Soil moisture (z0) are randomly distributed.

  • Soil moisture (z) can show aggreated patterns.

  • Ridge sites have lower soil moisture

  • Points indicate mean trait values for each grid (local community).
  • Ridge sites have greater trait values (e.g., thick dense leaves to grow well in dry conditions)
  • Is this because those traits are favored in ridge sites?
  • Is this just because neighbors have similar trait values (spatial autocorrelation)?

3 Models

\[ \left[ \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \end{array} \right] = \rho \left[ \begin{array}{rrrrrr} 0 & 1/3 & 0 & 1/3 & 1/3 & 0 \\ 1/3 & 0 &&&& \\ 0 && 0 &&& \\ 1/3 & & & 0 && \\ 1/3 & & & & 0 & \\ 0 & & & & & 0 \end{array} \right] \cdot \left[ \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \end{array} \right] + \epsilon \]

\[ y_1 = \rho * 1/3 (y_2 + y_4 + y_5) + \epsilon \]

In general, spatial autocorrelaiton models take this kind of form.

\[ Y = X \beta + \rho W Y + \epsilon \]

  • \(Y\): Dependent variable (\(N \times 1\), e.g., CWM)
  • \(X\): Independent variable (\(N \times k\), e.g., ridge (1) or valley (0))
  • \(\beta\): Model parameters (\(k \times l\))
  • \(\rho\): Spatial autocorrelaiton parameter
  • \(W\): Spatial weight matrix (\(N \times N\))
  • \(\epsilon\): errors ($N 1)

4 R examples

  • Let’s apply spatial autocorrelaiton models for the above violin plot example.
  • First, we need to prepare a spatial weight matrix (W).
  • Then, we apply some R functions to fit \(Y = X \beta + \rho W Y + \epsilon\)

4.1 Spatial weight matrix

  • Read the dataset.
d <- read_csv("./data/trait_env.csv")
d
## # A tibble: 625 x 4
##        x     y hab     trait
##    <dbl> <dbl> <chr>   <dbl>
##  1     1     1 ridge   2.93 
##  2     2     1 valley -0.819
##  3     3     1 valley -2.28 
##  4     4     1 ridge   4.99 
##  5     5     1 valley -0.996
##  6     6     1 valley -0.574
##  7     7     1 ridge   1.46 
##  8     8     1 ridge   0.596
##  9     9     1 ridge   1.26 
## 10    10     1 ridge   1.35 
## # … with 615 more rows
(nb <- cell2nb(length(unique(d$x)),length(unique(d$y)), type = "queen"))
## Neighbour list object:
## Number of regions: 625 
## Number of nonzero links: 4704 
## Percentage nonzero weights: 1.204224 
## Average number of links: 7.5264
(W <- nb2listw(nb, zero.policy = TRUE, style = "W"))
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 625 
## Number of nonzero links: 4704 
## Percentage nonzero weights: 1.204224 
## Average number of links: 7.5264 
## 
## Weights style: W 
## Weights constants summary:
##     n     nn  S0       S1       S2
## W 625 390625 625 169.8933 2510.443

4.2 Fit

  • Intercept only
fit1 <- spautolm(trait ~ 1, listw = W, data = d)

summary(fit1)
## 
## Call: spatialreg::spautolm(formula = formula, data = data, listw = listw, 
##     na.action = na.action, family = family, method = method, 
##     verbose = verbose, trs = trs, interval = interval, zero.policy = zero.policy, 
##     tol.solve = tol.solve, llprof = llprof, control = control)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -4.9026232 -0.4330284  0.0047709  0.4380002  5.1772142 
## 
## Coefficients: 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.024362   0.094447  0.2579   0.7965
## 
## Lambda: 0.60843 LR test value: 128.96 p-value: < 2.22e-16 
## Numerical Hessian standard error of lambda: 0.045821 
## 
## Log likelihood: -857.8677 
## ML residual variance (sigma squared): 0.85481, (sigma: 0.92456)
## Number of observations: 625 
## Number of parameters estimated: 3 
## AIC: 1721.7

  • lambda (rho) = 0.6084325 indicates a positve spatial autocorrelaiton

  • Model with an environmetal predictor

fit2 <- spautolm(trait ~ hab, listw = W, data = d)

summary(fit2)
## 
## Call: spatialreg::spautolm(formula = formula, data = data, listw = listw, 
##     na.action = na.action, family = family, method = method, 
##     verbose = verbose, trs = trs, interval = interval, zero.policy = zero.policy, 
##     tol.solve = tol.solve, llprof = llprof, control = control)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.30126 -0.41144  0.02866  0.40400  4.36665 
## 
## Coefficients: 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  0.702901   0.052053  13.504 < 2.2e-16
## habvalley   -1.370559   0.065657 -20.875 < 2.2e-16
## 
## Lambda: 0.25345 LR test value: 12.031 p-value: 0.00052329 
## Numerical Hessian standard error of lambda: 0.070405 
## 
## Log likelihood: -717.2709 
## ML residual variance (sigma squared): 0.57583, (sigma: 0.75883)
## Number of observations: 625 
## Number of parameters estimated: 4 
## AIC: 1442.5
  • lambda (rho) = 0.2534548 indicates a weak positve spatial autocorrelaiton
  • Negative habvalley indicates valley sites have negative effects on trait values compared to ridge sites even after controlling spatial autocorrelaiton.

4.3 Exercise

  • Can you make a raster figure for \(\epsilon\)?
    • fit1$fit$residuals is a vector of \(\epsilon\).
    • What does this figure mean?
  • Can you make a violin plot or box plot for \(\epsilon\)?
    • geom_violin, geom_boxplot
    • What does this figure mean?